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During last two decades it has been discovered that the statistical properties of a number of 
microscopically rather different random systems at the macroscopic level are described by the same 
universal probability distribution function which is called the Tracy-Widom (TW) distribution. 
Among these systems we find both purely methematical problems, such as the longest increasing 
subsequences in random permutations, and quite physical ones, such as directed polymers in random 
media or polynuclear crystal growth. In the extensive Introduction we discuss in simple terms these 
various random systems and explain what the universal TW function is. Next, concentrating on 
the example of one-dimensional directed polymers in random potential we give the main lines of 
the formal proof that fluctuations of their free energy are described the universal TW distribution. 
The second part of the review consist of detailed appendices which provide necessary self-contained 
mathematical background for the first part. 
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I. INTRODUCTION 

Everyone knows the Gaussian distribution function. Whenever we are dealing with a system containing independent 
random parameters, its macroscopic characteristics (according to the central limit theorem) are described by the 
Gaussian distribution. This kind of universal behavior is trivial, and not so much interesting. On the other hand, 
every non-trivial system usually requires individual consideration, and although there are lot of universal macroscopic 
properties among microscopically different systems (e.g. scaling and critical phenomena at the phase transitions) 
until very recently no one would expect to have a universal function (different from the Gaussian one) which would 
describe macroscopic statistical properties of a whole class of non-trivial random systems. 

Originally the solution of Tracy and Widom [l[ were devoted to rather specific mathematical problem, namely the 
distribution function of the largest eigenvalue of N x N Hermitian matrices (Gaussian Unitary Ensemble (GUE)) in 
the limit N — > oo. Nowadays we have got rather comprehensive list of various systems (both purely mathematical and 
physical) whose macroscopic statistical properties are described by the same universal Tracy- Widom (TW) distribution 
function. These systems are: the longest increasing subsequences (LIS) model [2\ (Section LA) zero-temperature lattice 
directed polymers with geometric disorder || the polynuclear growth (PNG) system (Section LB) the oriented 
digital boiling model the ballistic decomposition model @, the longest common subsequences (LCS) Q, the one- 
point distribution of the solutions of the KPZ equation [1] (which describes the motion of an interface separating two 
homogeneous bulk phases) in the long time limit 0, [nj, and finally finite temperature directed polymers in random 
potentials with short-range correlations [LlT - [l4| ]. It should be noted that directed polymers in a quenched random 
potential have been the subject of intense investigations during the past three decades (see e.g. (Isj). Diverse physical 
systems such as domain walls in magnetic films [lq ] , vortices in superconductors [l7| , wetting fronts on planar systems 
[18j . or Burgers turbulence [191 ] can be mapped to this model, which exhibits numerous non-trivial features deriving 
from the interplay between elasticity and disorder. 

The rest of this Introduction is devoted to the discussion of several random statistical systems with explanations 
in very simple terms of what the TW distribution describes in them. Namely, we will consider the combinatorial 
model of the longest increasing subsequences (section LA), the polynuclear crystal growth model (section LB), and 
one-dimensional directed polymers in random potential (section I.C). Besides in section I.D the main ides of the 
replica method (used in the present approach) will be described, and finally in section I.E the definition and the main 
properties of the TW distribution function will be given. 

Sections II and III are devoted to the exact solution of the one-dimensional directed polymers problem. In particular, 
in section II the main ideas of this solution as well as its methodological tools are described. Section III contains the 
main lines of the derivation of the TW distribution function for the free energy fluctuations in one-dimensional directed 
polymers with (5-correlated random potential. The second part of this review contains several technical appendices 
containing all necessary mathematical tools (used in the previous sections) which hopefully makes the whole paper 
to be self-contained. 



A. Combinatorics 

We start with purely mathematical "toy" model which, as we will see later, is closely related with physical problems 
of polynuclear crystal growth (section LB) and one-dimensional directed polymers (section I.C). This combinatorial 
problem of statistical properties of the longest increasing subsequences (LIS) was formulated long time ago by Ulam 
20], and hence it is often called the Ulam's problem. Let us consider a sequence of N integers {1,2,..., N}. Then, for 
an arbitrary permutation of these integers we have to find all possible increasing subsequences, and among them the 
length Zjv of the longest ones should be defined. As an example let us consider the case N = 11 and take a particular 
permutation 

{3,5,10,1,9,6,8,4,7,11,2}. (1.1) 

This permutation exhibits many different increasing subsequences (such as {3,5,10,11}, {1,9,11} etc.), and among 
them the longest ones are {3,5,6,7,11} and {3,5,6,8,11}. In other words, for this particular permutation In = 5. 
Simple graphical representation of this permutation problem is shown in Figure 1. Here the set of 11 bold dots inside 
the (12 x 12) square represents the permutation, cq. (|l.l[) : for every integer in the x direction one associates one and 
only one permuted integer in the y direction (for x = 1 one gets y = 3, for x = 2 one gets y — 5, etc.). All possible 
increasing subsequences of this particular permutation are obtained by drawing all possible directed paths connecting 
the origin (0, 0) with the right up corner (12, 12) of this square, which are passing over the internal bold dots. Directed 
path means that only "right-and-up" movements are allowed when going from one dot to another. For example, from 
the point (6, 6) one can jump only to the points (7, 8), (9, 7) and (10, 11). In other words, when going from the origin 
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to the right up corner (12, 12) both x and y coordinates can only increase at every step. In terms of these rule, the 
longest increasing subsequence is described by the directed path which goes over the maximum possible number of 
dots. Note that for a given permutation the longest increasing subsequence is not necessary unique. In the considered 
example in addition to the subsequence {3, 5, 6, 7, 11} shown in Figure 1 by the dotted line, there exist another one 
namely {3, 5, 6, 8, 11}. 
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FIG. 1: Geometrical representation of the permutation, eq, (|l.l[) , for N = 11. The dotted line corresponds the longest increasing 
subsequence {3, 5, 6, 7, 11} 



Considering all N\ possible permutations with equal probability one find that the length In varies from permutation 
to permutation being the random variable. The question is what are the statistical properties of the random quantity 
In- First it has been shown that in the limit of large N, the average value of the longest increasing subsequence, In, is 
proportional to vN, namely In ~ cyN, where the constant c = 2 21]. Moreover, in the limit N — » oo, the quantity 
In in order \/~N is self averaging: limjv->oo In/\/N = 2 (in other words, the distribution function of the ratio In/VN 
in the limit N — > oo shrinks to the S- function). 

Note that at the qualitative level one can easily understand why the typical value of In must be proportional to \/N. 
Indeed, for large N a generic permutation of the numbers {1,2,..., N} in terms of the permutation matrix of Figure 
1 will be represented by a uniform distribution of N dots inside the N x N square. Thus the density of the dots will 
be proportional to 1/N while the typical distance between neighboring dots must be proportional to y/~N. It means 
that the typical number of dots on a diagonal-like path of the length ~ N must be proportional to N/ ViV = y/N. 

However, the main interest in this system is not the typical value of the longest increasing subsequences, but their 
fluctuations. Recently it has been shown 0, [22| that in the limit of large N the fluctuations of In scale as TV 1 / 6 , 
namely, 

l N 2\fN + N^ 6 s (1.2) 

where the random quantity s is described by the universal N- independent distribution function Ptw(s), which is the 
Tracy- Widom distribution (see section I.E): 

J i ? 00 Pr0b ( k ^ = S ) = Ptw{3) (L3) 

It turns out that the above purely mathematical "toy" model is equivalent to the physical model of (2+1)- 
dimensional polynuclear crystal growth, where the TW distribution describes the fluctuations of the number of the 
crystal mono-layers (see next subsection). 
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B. Polynuclear crystal growth 



It turns out that the mathematical " toy" model considered above is equivalent to the physical model of the crystal 
growth with randomly located nucleation centers. This is the model of polynuclear crystal growth (PNG) which 
describes the growth of the two-dimensional crystal monolayers in (2+1) dimensions. 

Let us consider again Figure 1 where the bold dots will be assumed to represent the nucleation centers. The crystal 
layers growth takes place in the vertical direction (toward the reader) according to the following rules. From each 
nucleation center we draw the monolayer level step straight line in the horizontal direction to the right and in the 
vertical direction up, until these lines meet with the other lines starting from the other centers. In this way we are 
getting the monolayer "terraces" which mount from the left-down to the right-up corner of the square, (see Figure 2). 
One can easily see that for a given random positions of the nucleation centers (for a given permutation in the previous 
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FIG. 2: (2+l)-dimensional crystal terraces with the "nucleation centers" corresponding to the permutation model of Figure 1 



LIS problem) the number of terraces /ijv is just equal to the longest increasing permutation In in the combinatorial 
problem considered in the previous subsection. For a given value of N, depending on the actual configurations of 
the nucleation centers inside the N x N square, the number of the monolayer terraces is the random quantity, and 
in the limit N —> oo its statistics is described by the TW distribution, eq. (jl.3[> [4J . At present the study of various 
modifications of PNG model formulated above is the vast field of research (see e.g. [23j]). It is also interesting to note 
that quite recently the existence of TW distribution in the PNG-like systems has been confirmed experimentally [2~3 | . 



C. Directed polymers 



It is clear that when the size of the square N is large the presence of the background lattice is not essential. In this 
case instead of considering the problem in terms of permutations one can introduce a homogeneous distribution of 
dots inside continuous square (Figure 3). Let us introduce here the diagonal "time" axes which goes from the origin 
(left down corner) to the right up corner of the square. The directed polymer here is the line which starts at the 
origin and arrives to the right up corner "jumping" over the dots in such a way that its time coordinates increases at 
every jump. For a given (random) configuration of dots among many possible directed polymer trajectories we have 
to choose the ones which contains the maximum number of dots. In this formulation the problem of statistics of the 
length of such directed polymers looks somewhat different from the LIS problem considered above. It can be shown 



FIG. 3: Directed path through randomly distributed dots 



however, that in the limit of large times t and large N these two problems become equivalent [2]. For a given (fixed) 
density p of dots instead of the total number of dots N one can measure the length l(t) of the polymer in terms of 
the size of the square t. Since p = 2N/t 2 , assuming that p is a parameter which is of the order of one, we note that 
t oc \/N. In this case we find that the fluctuations of the length of the polymers scales as i 1 / 3 , and instead of eq. (11.3[) 
we get 

In other words in the thermodynamic limit, t — > oo all the systems considered above appear to be equivalent to each 
other which doesn't look so much surprising if we compare Figures 1-3. 

In statistical physics one defines random directed polymers in somewhat different way. Let us introduce a square 
lattice in which discrete "time" t=l,2,...,Lis now horizontal (Figure 4). The vertical direction is described by the 
discrete parameter <f> = 0,±1,±2, ...,±M. At every lattice site (<j>,t) we place random quantities (random potential) 
V((f>, t) and assume that they are described by independent Gaussian distributions: 

<t>,t 

The parameter u defines the typical strength of the random potentials V{4>, t), which according to eq. (jl.5l) arc 
uncorrelated 



V{<t>,t)V{<l/,t') = u5^5 t ,t> (1.6) 

and have zero mean value, V(<p,t) = (the horizontal line denotes the averaging with the distribution, eq. (|1.5|l ). 

The directed polymer here is the path which starts at the origin and goes over the lattice sites to the right end of 
the system. At every time step, t — > t + 1 the polymer trajectory <j>{t), can deviate up or down by one step or may not 
deviate at all: 4>(t + 1) = <p(t) + a(t + 1), where a = ±1, 0. Assuming that the polymer is a kind of elastic string we 

can introduce "elastic" (positive) energy oc \<j>(t + 1) — 4>(t)~\ for every polymer's deviation. In this way for a given 
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FIG. 4: Directed polymer on a square lattice 



trajectory <p(t) we can associate the following energy: 



L-l 



This expression contains two competing contributions: the first (elastic) terms are trying to make the trajectory as 
horizontal as possible, while the second ones are forcing the trajectory to deviate in the search for the most negative 
values of the random potentials. For a given (random) configuration of the potentials V(4>, t) the optimal trajectory 
4>*(t) is defined by the minimum of the total energy, eq. (|1.7j) . 



E(L) = mm j £ [- [cf>(t + 1) - <f>(t)] + V(tf>(t),t)\ j (1.8) 

Being the function of the random potential V(cj>,t) this quantity is also random, and it could be considered as the 
distant analog of the (random) length of the directed polymers in the previous example shown in Figure 3. An essential 
difference is that in this latter case the elastic terms are absent, while the (negative) contribution of the potential 
energy is associated with a fixed energy Vq carried by the dots which (unlike the Gaussian random potentials V{<j>, t)) 
are geometrically random. One more important difference is that unlike the directed polymer in Figure 4, which is 
defined by the local in space one-step wandering, the trajectory in Figure 3 can jump any distance from dot to dot 
(not necessary between neighboring dots). Thus, a priori there are no reasons to expect that these two quantities, 
E(L), ea. dl.8p . and l(t) in the example of Figure 3 would have the same statistical properties. If, nevertheless, we 
would suppose that at least in the thermodynamic limit, t ~ L — > oo, these two types of systems become equivalent 
we would have to expect that E(L) ~ /o L + L 1 / 3 s, where /o is the bulk (selfaveraging) energy density and the 
random parameter s is described by L-independent TW distribution function. 

The system defined by the Hamiltonian (jl.7l) is the usual one-dimensional statistical system containing quenched 
disorder. The fact that the leading contribution to its ground state energy E(L) is proportional to the system size L 
can be explained in very simple way. Indeed, in the first approximation, to minimize the energy at every time step 
among three possible options ("up" /'horizontal" and "down") the trajectory of Figure 4 can choose the site where the 
value of the random potential V(tfi, t) is lower (in this approximation we neglect the presence of non-local phenomena 
when the optimal trajectory chooses locally unfavorable option to gain globally more favorable energy). In this way 
the second term in eq. (|1.7p provides the contribution which is proportional to —y/uL (and not \J uL, as it would be 
for an arbitrary trajectory 4>(t)). Since the contribution of the first (elastic) term in eq. (|1.7|) is also proportional to L, 
we find that in the leading order in L, the energy of the optimal trajectory E(L) ~ —{const) L, and moreover, we 
can be sure that this contribution is negative since the energy of the optimal trajectory in the absence on the random 
potentials is zero (it is just the straight horizontal line), while the presence of the random potential can only lower 
the energy. On the other hand the fact that finite-size corrections in this system are of order L 1 / 3 (and not of order 
L 1 / 2 , as one could naively expect) is highly non-trivial phenomenon which is very difficult to explain in simple terms. 

The lattice model as it is introduced above, eas. (ll.7[) - (|1.8|i . is essentially the zero-temperature system, as we are 
dealing here with the optimal (global minimum) trajectories only. It is clear that the search for the global minimum 
configurations in eq. (|1.8[) is highly non-trivial task, as it can not be done via the local in time algorithms. On the 
other hand, as it often happens, one can make life much easier if one consider more general (i.e. more complicated) 
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problem. Namely, let us introduce finite temperature T in the system so that in addition to the quenched disorder 
fluctuations we would have the contributions of the thermal fluctuations produced by the trajectories away from the 
global minima ones. In terms of this generalization instead of the global minimum energy E(L), eq. (|1.8[) . we would 
get the free energy: 



F(L,T) = -Tln[£exp(-iflMt)] 

0(t) 



(1.9) 



where the expression under the logarithm is the partition function in which the summation goes over all trajectories 
starting at the origin. In the case the global minimum trajectory is unique (which is usually the case in the large 
system) the energy defined in eq. (|1.8|) is obtained by taking the zero temperature limit: E(L) = lirriT^o F(L, T). The 
finite temperature lattice model defined by eqs. (ll.7[ ). (|1.9[) is sufficiently simple for numerical investigations (see e.g. 
25]) which in particular allows to demonstrate the existence of the free energy fluctuations scaling ~ L 1 / 3 . On the 
other hand, the problem defined on a lattice is very hard for analytical studies. Since the phenomena of interest are 
taking place at large system sizes, one may hope that it would be sufficient to consider the system in the continuous 
limit where the presence of a lattice become irrelevant. Continuous limit generalization of the Hamiltonian, eq. (ll.7[) 
is straightforward. Assuming that the lattice spacing goes to zero and changing the finite differences in its first term 
by the gradients we get 





r#(r)l 




dr 



v{4>{t),- 



(1.10) 



where, as before, the disorder potential V(4>, r) is Gaussian distributed and uncorrelated. Instead of eq. (|1.6[) in the 
continuous limit we get 



V((t>,T)V(<t)',T') = u8{<P~(j)')5{T-T') 



The partition function of this system is now defined in terms of the functional integral: 

/+oo p<j>(L)=4> L 
d<p L / V[4>(t)] 



(1.11) 



(1.12) 



where f3 — 1/T is the inverse temperature and the integration goes over all trajectories starting at the origin (t = 0) 
and having free boundary conditions at r = L. In this way instead of the lattice trajectory of Figure 4 we are getting 
a continuous trajectory shown in Figure 5. Although, as we will see later, the continuous model defined above is ill 
defined at short distances, as far as its long-time behavior is concerned it is much better treatable analytically. 




FIG. 5: Continuous elastic string in a random potential, eq. (|1.10[ ) 



First of all we can note that in the absence of the random potential in the Hamiltonian, ea. (|1.10|i . the system 
describes simple thermal diffusion. Indeed, the probability that at time r = L the trajectory arrives to the point 
<f>(L) = cj)L is given by the partition function 







<f>(L)=<p L 



<Ko)=o 



(1.13) 
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Simple Gaussian integration (with the proper choice of the integration measure of the functional integral) yields 

z ° fe) = 7lb ""fir} (L14) 

In other words the typical deviation < 4>l > of the trajectory due to the thermal fluctuations growth as VT L 1 / 2 
(and it shrinks to zero at T = 0). On the other hand, in the presence of the random potential the two terms of the 
Hamiltonian (|1.10[> must balance each other. For a given value of the typical deviation the contribution of the 
elastic term can be estimated as <fijL. Thus, if in the presence of disorder the free energy fluctuations of this system 
are scaling as L 1 / 3 (see e.g. |26l - |29| ) we can conclude that the typical value of the trajectory deviations due to the 
action of the random potentials must grow as (f>L ~ L 2 / 3 which is much faster than the pure thermal diffusion scaling 

In what follows we are going to study the system defined by eqs. (|1.10p . (|l.ll[) in all detail. It turns out that regardless 
of essential differences in the definition of this system and the previous models discussed above in sections l.A-l.C, 
in the thermodynamic limit all these models becomes equivalent. The central result which will be proved in the next 
sections is in the following. In the limit L — > oo the free energy of the considered system can be represented as 

F = f L + cL 1 ' 3 f (1.15) 

where /o is the linear free energy density, the constant c = ^(^m 2 ) 1 / 3 , and the random quantity / ~ 1 (just like the 
quantity s in eqs. (|1.2[) - (|1.4l0 is described by the universal TW distribution function (see section I.E) 



D. Replica method 

The replica method is widely used in the studies of systems containing quenched disorder (see e.g. [3(1 HH). For 
simplicity let us consider the string </>(t) with the zero boundary conditions: 0(0) = 4>{L) = 0. The partition function 
of a given sample described by the Hamiltonian, eq. (11.101) , is 

r 4>{L)=o 

z[v] = / v[4>{t)\ e-^< y i (i.i6) 

J0(O)=O 

On the other hand, the partition function is related to the total free energy F[V] via 

Z[V] = exp(-/3F[V]) (1.17) 

The free energy F[V] is defined for a specific realization of the random potential V and thus represent a random 
variable. Taking the iV-th power of both sides of this relation and performing the averaging over the random potential 
V we obtain 



Z N [V] = Z[N, L] = exp(-/3NF[V]) (1.18) 

where the quantity in the lhs of the above equation is called the replica partition function. Substituting here the free 
energy in the form F = JqL + cL 1 / 3 /, cq. (|1.15p . and redefining the partition function 

Z[N, L] = Z[N, L] e~ pNfoL (1.19) 

we get 



Z[N, L] = exp(-A7V/) (1.20) 

where A = ficL 1 / 3 . The averaging in the rhs of the above equation can be represented in terms of the distribution 
function Pl{J) (which depends on the system size L). In this way we arrive to the following general relation between 
the replica partition function Z[N,L] and the distribution function of the free energy fluctuations Pl(/): 

r+oo 

Z[N,L] = / dfP L (f)e- XN f (1.21) 



Of course, the most interesting object is the thermodynamic limit distribution function P*(/) = limi^oo Pl(I) which 
is expected to be the universal quantity. The above equation is the bilateral Laplace transform of the function 
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Pl(/), and at least formally it allows to restore this function via inverse Laplace transform of the replica partition 
function Z[N, L]. In order to do so one has to compute Z[N, L] for an arbitrary integer N and then perform analytical 
continuation of this function from integer to arbitrary complex values of N. This is the standard strategy of the replica 
method in disordered systems where it is well known that very often the procedure of such analytic continuation turns 
out to be rather controversial point [32l [33j . Even in rare cases when the derivation of the replica partition function 
Z(N) — Z N can be done exactly, its further analytic continuation to non-integer N appears to be ambiguous. The 
classical example of this situation is provided by the Derrida's Random Energy Model in which the momenta Z(N) 
growths as fast as exp(iV 2 ) at large N, and in this case there are many different distributions yielding the same values 
of Z(N), but providing different values for the average free energy of the system [34[. In our present system the 
situation is even worse because, as we will see later, the replica partition function growth here as exp(N 3 ) at large TV, 
and in this situation its analytic continuation from integer to non integer N would be rather problematic point. It 
turns out, however, that in our present case it is possible to bypass the problem of the analytic continuation if instead 
of the distribution function P*(f) one would study its integral representation 



W(x) = / dfP^f) (1.22) 

J X 

which gives the probability to find the fluctuation / bigger that a given value x. Formally the function W{x) can be 
defined as follows: 

W(x) = lim V V ' exp(AiVir) Z N (1.23) 

N=Q 

00 (—l) N 

= , lim E S^r- e M*Nx - \Nf) 

A— >oo * — ' iV ! 

N=0 



= lim exp[-exp(A(x-/))] = 9(f - x) 

A— >oo 

Thus, the probability function W{x) can be computed in terms of the above replica partition function Z\N, L] by 
summing over all replica integers 

W(x) = lim V 1 ' exp(XNx) Z[N, L] (1.24) 

Of course, keeping in mind that Z[N,L] ~ exp(iV 3 ) at large N, we see that the above series is not that innocent. 
Here in accordance with the troubles conservation law instead of the problem of analytic continuation we are facing 
formally divergent series. Nevertheless, it can be shown that this sign alternating series can be regularized in the 
standard way (similarly to formally divergent sign alternating series X]fcLo( — l) fc a fc = (1 + a)^ 1 -* which at \a\ > 1 is 
defined as the analytic continuation from the region \a\ < 1). This eventually allows to prove that the thermodynamic 
limit function W(x), eq. (|1.24j) . is defined by the universal Tracy- Widom distribution function. 



E. Tracy- Widom distribution function 



Originally the Tracy- Widom distribution function has been derived in the context of the statistical properties of 
the Gaussian Unitary Ensemble (GUE) of random Hermitian matrices [1]. GUE is the set of N x N random complex 
Hermitian matrices (such that = G* { ) whose elements are drawn independently from the Gaussian distribution 



P[G] = B N exp{-irr(G 2 )} 



(1.25) 



where Bn is the normalization constant. The joint probability density of N eigenvalues {Ai, A2, Xn} of such 
matrices has rather compact form 



V[Xi,X 2 ,...,X N ] 



N 



exp 



N 

{-!>•} 



(1.26) 



where Cn is the normalization constant. Using this joint probability density one can calculate various averaged 
characteristics of the eigenvalue statistics. For example, one can introduce the average density of the eigenvalues 
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p(\,N) = ^5^i_i(5(A — Aj)) where the averaging ( ) is performed with the probability distribution, eq. (|1.26l) . 

Using the symmetry of this distribution one gets 



p(X,N) = 

It can be shown [35[ that in the limit of large N 

P (X,N) 



N 

n 

i=2 



+00 



V[X,X 2 ,...,X N ] 



(1.27) 



Nn 2 



A2_ 

27V 



(1.28) 



We see that on average the eigenvalues lie within the finite interval [— y/2N < X < y/2N\ where, according to eq. (|1.28|l . 
their density has the semi-circular form. This is one of the central results of the random matrix theory which is called 
the Wigner semi-circular law. In particular this result tells that on average the maximum eigenvalue X max is equal 
to \/2N . However, at large but finite N the value of X max is the random quantity which fluctuates from sample to 
sample. One may ask, what is the full probability distribution of the largest eigenvalue X max 7 This distribution can 
be computed in terms of the general probability density, eq. (ll.26[) . Introducing standard notations of the random 
matrix theory we define the function F2(s) = Prob[A ma:E < si which gives the probability that X max is less than 
a given value s (in these notations the functions -Fi(s), F 2 (s) and F^s) denote the probability distributions of the 
largest eigenvalues in the Gaussian Orthogonal Ensemble (GOE), Gaussian Unitary Ensemble (GUE) and Gaussian 
Symplectic Ensemble (GSE) correspondingly [36]). By definition 



F 2 (s) = 



N r s l r s 

II / ^ ^^2,...^]= / dXPrwW 

A — A J —OO J —OO 



(1.29) 



It is this problem which has been solved by Tracy and Widom in 1994 1]. It has been shown that at large N the 
typical fluctuations of X max around its mean value y/2N scale as TV -1 / 6 , namely (c.f. eq. (|1.2l0 



A r , 



V2N 1 / 6 



(1.30) 



where the random quantity s is described by iV-independent distribution Ptw(s) = dF 2 {s)/ds. The function F 2 (x), 
has the following explicit form 



F 2 (s) = exp^-^ 



dt (t-s) q 2 {t) 



(1.31) 



Ptw(s) = -i-F 2 (s) = exp 
ds 



dt (t- s)q 2 (t) 



dtq 2 (t) 



where the function q(t) is the solution of the Panleve II equation^ 1 ), 

q" = tq + 2q 3 



(1.32) 



(1.33) 



with the boundary condition, q(t — > +00) ~ Ai(t). The shape of the function Ptw{s) is shown in Figure 6. Note 
that the asymptotic tails of this function are strongly asymmetric. While its right tail coincides the Airy function 
asymptotic Ptw(s — > +00) ~ exp[— |s 3 / 2 ] , the left tail exhibits much faster decay Ptw( s —00) ~ exp[— ^|s| 3 ] 



[1] There exist six Panleve differential equations which were discovered about a hundred years ago[37l] (for the recent review see e.g. [38l|). 
It is proved that the general solutions of the Panlevee equations are transcendental in a sense that they can not be expressed in terms of 
any of the previously known function including all classical special functions. At present the Panleve equations have many applications 
in various parts of modern physics including statistical mechanics, plasma physics, nonlinear waves, quantum field theory and general 
relativity 
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^6 -4 -2 "0 2 

FIG. 6: Tracy- Widom distribution function Ptw(x), eg. (|1.32l) . 



II. MAPPING TO QUANTUM BOSONS 



Explicitly, the replica partition function, Eqs. (|1.18|) . f| 1 . 16[) . of the system described by the Hamiltonian, Eq. (|1.10j) . 

is 



N 



4>a(L)=0 



Z(N,L) = U / 2?0 a (r)exp 



a=1 J0 a (O)=O 



L N 



a=l 



(2.1) 



Since the random potential V[4>, r] has the Gaussian distribution the disorder averaging (...) in the above equation is 
very simple: 



exp 



r L n r 2 r N 

-p / drY,V[Mr),r] = exp tL / / drdT ' £ V[<p a (r), r]V[<Mr'), A 

J ° a=l J L J J ° a,b=l 



(2.2) 



Using Eq. (jl.lip we find: 



Z(N,L) = U 



0o(L)=O 



a=1 J4> a (o)=o 



2?</> Q (r) exp 



JV N 



a=l 



„b=l 



(2.3) 



It should be noted that the second term in the exponential of the above equation contain formally divergent contri- 
butions proportional to 6(0) (due to the terms with a = b). In fact, this is just an indication that the continuous 
model, Eqs. (|1.10p - (|l.lip is ill defined as short distances and requires proper lattice regularization. Of course, the 
corresponding lattice model Eqs. (|1.6p - (|1.7p contains no divergences, and the terms with a = b in the exponential of 
the corresponding replica partition function would produce irrelevant constant ^L(3 2 uN5(0) (where the lattice version 
of 6(0) has a finite value). Since the lattice regularization has no impact on the continuous long distance properties 
of the considered system this term will be just omitted in our further study. 
Introducing the iV-component scalar field replica Hamiltonian 



1 f L ( N N 
Hn[cI>} = - / dr[J2[drMr)] 2 ~ PuY,S[Mr) 

1 J ° \a=l a^b 




(2.4) 



for the replica partition function, Eq. 



we obtain the standard expression 

N r a (L)=O 



Z(N,L) = Yl 



a=1 J<t>*(o)=o 



(2.5) 



where = {<f>%, . . . , ((>n}- According to the above definition this partition function describes the statistics of N 
6- interacting (attracting) trajectories </> a ( T ) all starting (at r = 0) and ending (at r = L) at zero: <j> a (0) = 4> a (L) = 
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In order to map the problem to one-dimensional quantum bosons, instead of the above replica partition function, 
cq. (j2.5[) , let us introduce more general object 

*(x;i) = TT/ V<p a (r) e-? H »M (2.6) 

which describes N trajectories (j> a (r) all starting at zero (4> a {0) — 0), but ending at r = t in arbitrary given points 
{xi, ...,xn}. One can easily show that instead of using the path integral, ^(xji) may be obtained as the solution of 
the linear differential equation 

N N 

Wx;t) = 2flE^ a *(x;*) + -P 2 uJ2S(xa~x b )^(x;t) (2.7) 

" 0=1 a^b 

with the initial condition 

*(x;0)=n^ =1 <y(a: a ) (2.8) 
One can easily see that Eq. (|2.7p is the imaginary-time Schrodinger equation 

-ft*(x;t) = H*(x;t) (2.9) 

with the Hamiltonian 

JV JV 
P o=l 

which describes iV bose-particles of mass /3 interacting via the attractive two-body potential —j3 2 u5(x). The original 
replica partition function, Eq. (|2.5[) . then is obtained via a particular choice of the final-point coordinates, 

Z(N,L) = *(0;X). (2.11) 

Historically the main interest in the studies of such type of systems was devoted to the quantum bosons with 
repulsion. It is for the case of repulsive interactions the free energy of the N particle system reveals "correct" 
extensive behavior, F oc N. The eigenfunction of iV-particle Hamiltonian Eq. (|2.10[) . with repulsive interactions, 
u < 0, have been derived by Lieb and Liniger in 1963 [39[ (for details see Appendix A, as well as Refs. [I(||4lj]). The 
system of attractive bosons remained much less studied. The free energy of such system reveals "bad" thermodynamic 
limit behavior, namely F cx —N 3 . Besides, as will be shown later, the structure of the eigenstates of such system is 
much more complicated compared to the case of repulsion. The spectrum and some properties of the eigenfunctions 
for attractive (u> 0) one-dimensional quantum boson system have been derived by McGuire [42[ and by Yang [43[ 
(see also Ref. @, HU). Detailed structure and the properties of these wave functions are described in Appendix B. 

A generic eigenstate of such system consists of M (1 < M < N) "clusters" {f2 Q } (a — 1, M) of bound particles. 
Each cluster is characterized by the momentum q a of its center of mass motion, and by the number n a of particles 
contained in it (such that n a — N). Correspondingly, the eigenfunction ^q,n (%x, %n) of such state is 

characterized by M continuous parameters q = (q±, <?m) and M integer parameters n = (ni, ...,tim) (see Appendix 
B2, Eq. (|B.27|) ). The energy spectrum of this state is 

1 M K 2 M 

" a=l " a=l 

where 

k = (3 3 u (2.13) 



A general time dependent solution ^(x, t) of the Schrodinger equation (|2.7p with the initial conditions, Eq 
be represented in the form of the linear combination of the eigenfunctions ^q^n (x) (Appendix B4, Eq. (IB.50p ): 

*(x,t) = E E' l' V * exp[-E M (q,n)t] (2.14) 



M=l n 
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where the summations over the integer parameters n a and the integrations over the momenta q a are performed in a 
restricted subspace, Eqs. (|B.42p - (|B.45p and (|B.51[) . which reflects the specific symmetry properties of the cigcnfunctions 

^q^n (x). Correspondingly, according to Eq. ()2.11j) . for the replica partition function of the original directed polymer 



problem one gets (Appendix B4, Eq. (fB759|) ) 



Z(N,L) 



N r M r+OQ , oo -, , M 

e in/ ^£ «fe 



M=\ 



N 



1*^(0) 



-E M (q,n)£ 



(2.15) 



where due to the symmetry of the function /(q, n) = l^q^CO)! 2 exp[— Em{<\, n) L\ with respect to permutations of 
all its M pairs of arguments (q a ,n a ) the integrations over M momenta q a can be extended to the whole space Rm 



while the summations over n a 's are bounded by the only constrain X)a=i n a = N (for simplicity, due to the presence 
of the Kronecker symbol S(J2 a n a , /V), the summations over n Q 's are extended to infinity). 



III. SOLUTION OF THE ONE-DIMENSIONAL DIRECTED POLYMERS PROBLEM 



Using the explicit form of the wave functions Vf^n (x), Eq. (|B.27j) . the expression in Eq. (|2.15j) for the replica partition 
function can be reduced to (Appendix B4, Eq. ([BT6T|) - (|B~62l) ) 

Z(N,L) =e~ pNLfo Z(N,X) 



(3.1) 



where /o = ^/3 4 u 2 — ln(/3 3 u) is the linear (selfaveraging) free energy density (cf. Eq. ([1.19[) ). and 



Z(N.L) = Nl 



+ OC 



dq 



exp 



, 2ttkN 

N r M oo r+00 



L , T 9 k 2 L , 
Nq 2 H N 

2(3 q 24,3 



M! 

M=2 L a=lnQ = 1 
A/ 



dg 



2irKn c 



/ M s A/ I 

\*=1 ' a<fl ?Q 



q a -qp- _ n /a)l 



A/ 



x exp 



.Zj ^ — > 2 Ki L ^ — *> -i 



2/3 



24/3 ^ 



(3.2) 



The first term in the above expression is the contribution of the ground state (M = 1), while the next terms (M > 2) 
are the contributions of the rest of the energy spectrum. 

The terms cubic in n a in the exponential of Eq. (|3.2|) can be linearised with the help of Airy function, using the 
standard relation (see Appendix C) 

r + oo 



exp 



dy Ai(y) exp(An) 



1 /3 

Redefining the momenta, q a = m/s/LJ p a and introducing a new parameter 



X(L) 



1 (L 



2\(3 



1/3 



1/3 



(3.3) 



(3.4) 



we get 

Z(N,X) = Nl 



N 



M=2 ' L a=l n a =l 



dydp 

4:TT\N 

r M oo 



Ai(y) e XN{y - p ) + 
+oc dy a dp a 



4TrXn a 



Ai{y a ) e 



M 

n 



|A(n Q - n^) - i(p a - pp)\ 



a <p \K n a +np) - i(p a -pp)\ 



(3.5) 



A/ 



After shifting the Airy function parameters of integration y a — > y Q +p„ the expression for Z(iV, A) becomes sufficiently 
compact: 

r+oc dydp 



Z{N,X) = Nl 



N 



M=2 L a =ln a =l 



47tAA^ 

r If OO 



Ai(y+p 2 )e XNy + 
+oc dy a dp a 



(3.6) 



47rAnQ, 



Ai(y a +P^) e 



2 \ pilaSa 



A/ 

n 



\X(n a - up) - i(p a -pp)\ 



a</3 \M n a +np) - i{p a -Pf})\ 



A/ 



8 [J2n a , N 
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Now, using the Cauchy double alternant identity 



U.a<a( a <* - ap)(b a ~ bp) 



n a ,s=i( a « -h 



= (_l)«(M-l)/2 det 



a a - bt 



= 1,...M 



(3.7) 



and introducing a a — p a — iXn a and b a = p a +iXn a , the product term in eq. (|3.6p can be represented in the determinant 
form: 



M 



n 



\\(n a - np) - i(p a -pp)\ 



p \X(n a + np) - i(p a -pp)\ 



M 



a=l 



det 



1 



Xn a — ip a + Xnp + ipp 



a,P=l,...M 



(3.8) 



Substituting now the expression for the replica partition function Z{N 1 X) into the definition of the probability function, 
eq. (|1.24[) . we can perform summation over N (which would lift the constraint X)q=i n a = -^0 ano - obtain: 



W(x) = Urn l+£ 



(-1) 



M 



M=l 



M\ 



M 

n 



dy a dp a 2 



2tt 



Mv.+pD £ (-l)—^""^^ 



det 



-An Q — ip Q + Xnp + ipp. 

(3.9) 

The above expression in nothing else but the expansion of the Fredholm determinant det(l — K) (see e.g. (46| . 
Appendix D) with the kernel 



k = K[(n,p);(n\p r )] = 



+ 00 



dyki{y + p 2 ){-l) n ^e 



n-lpAn(y+:r) 



1 



Art — ip + Xn' + ip' 



Using the exponential representation of this determinant we get 

\ 1 

W(x) — lim exp — > — 

M=l 



TtK 



M 



where 



TrK 



M 



' M 

n 

,a=l 



+ oo 



dy a dp a 
2tt 



Ai(y a + P 2 a ) x;(-ir- -i e 



-lnAn„(i/ a +i) 



n a = l 
1 



(Am - ipi + Xn 2 + ip2){Xn 2 - ip 2 + An 3 + ipz)...{Xn M - ipM + Xn% + ipi) 



Substituting here 



Xn a - ip a + Xn a+ i + ip a+ i 
one can easily perform the summation over n a 's. Taking into account that 

oo qXz 

lim r— 

a^oo i + e Xz 



du a exp[— (Xn a — ip a + Xn a+1 + ip a+1 )oj c 



lim V(-l)»- 1 e An * 

, — inn ' ^ 



A— J-OO 



we get 



M 



lim Tr K M = TT 

X^^i J- J- 



dy a dp a 
2ir 



duj a Ai(y a + p 2 a )0{y a + x - w a - w -i)G 



ip a (uJ a - 1 ) 



(3.10) 



(3-11) 



(3.12) 



(3.13) 



(3.14) 



(3.15) 



where by definition ujq = ljm- Shifting the integration parameters, y a — > y a — x + uj a + u tt -i and w Q — > Wq, + x/2, we 
obtain 



/.oo Z' + OO J />oo 

lim Tr K M = T[ dy a ^ doj a Ai(y a + p 2 a +uj a + Ua - ± ) & p ^ u 



-i) 



(3.16) 
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Using the Airy function integral representation, and taking into account that it satisfies the differential equation, 
Ai"(t) = t Ai(t), one can easily perform the following integrations (see Appendix C): 



+00 



dy I ^M{ V +p 2 +uj + uj')e ip{u - u ' ) = 2" 1 / 3 / dyAi(2 1 / 3 w + y)Ai(2 1 /V + 2; ) (3.17) 



J -00 27T Jq 



Ai(2 1 /3 W ) Ai'(2 1 /3 W ') _ Ai'(2 1 /3 W ) Ai(2 1 /3 W ') 



Redefining to a — > uj a 2 */ 3 we find 



A->oo 



lim TrK M = j ... j du l du 2 ...d^ M K A {uj 1 ,uj2)K A {uj 2 ,uj z )...K A {uj M ,uj 1 ) (3.18) 



where 



-x/2 2 /3 



K AM = ^Ml^iW (3 . 19) 



is the so called Airy kernel. This proves that in the thermodynamic limit, L — ¥ 00, the probability function W{x), 
cq. (|1.22|) . is defined by the Fredholm determinant, 

W{x) = det[l - iti] = F 2 (-a;/2 2/3 ) (3.20) 

where K A is the integral operator on [— x/2 2 / 3 ,oo) with the Airy kernel, eq. (13.19[) . This is the Tracy- Widom distri- 
bution function which has the following explicit form (see Appendix D): 

F 2 {s) = gx p(- dt (t- s) q 2 (t)^j (3.21) 

where the function q(t) is the solution of the Panleve II equation, q" = tq + 2q 3 with the boundary condition, 
q(t -> +00) - Ai(t). Note that according to Eqs.flTSl]), (|Q2]) and (13^01) . P*(x) = 2~ 2 / 3 P TW {-2- 2 / 3 x) 

IV. CONCLUSIONS 

The first breakthrough in the studies of one-dimensional directed polymers in random potential was due to the 
work of Kardar [29j , in which the problem has been reduced to A-particle system of quantum bosons with attractive 
interactions. By that time the very first idea was that in the thermodynamic limit it would be sufficient to take into 
account only the contribution of the ground state whose energy was well known. Indeed, for any integer N > 1 the 
contribution of the excited states in the limit L —¥ 00 are exponentially small compared to that of the ground state. 
In the framework of this approach it has been demonstrated that the free energy fluctuations grow as L 1 ! 3 while the 
typical value of the polymers deviations scale as L 2 ! 3 which, in particular, was in perfect agreement with numerical 
studies. 

However, more detailed investigations demonstrated that the above approach reveals serious pathologies. In par- 
ticular it turned out that the second cummulant of the free energy (F 2 — F ) appears to be identically equal to 
zero! This is possible only in two cases: either the quantity F is not random (which contradicts to the fact that its 
fluctuations scale as L 1 / 3 ), or the distribution function of the this quantity is not positively defined (which, of course, 
makes no physical sense) . Simple mathematical analysis demonstrated that the origin of this pathology is hidden in 
the replicas " magic operations" : on one hand, all the calculations are performed assuming that the replica parameter 
N (number of particles) is an integer N > 1, while on the other hand, in the thermodynamic limit L — > 00 the relevant 
values of the parameter N which defines the physical properties of the original random system appears to be in the 
region N — > 0. In other words, the replica method assumes analytic continuation of the result obtained for arbitrary 
integers N to the region N — ¥ 0. The problem is that, first, such analytic continuation is not always unambiguous 
(see e.g. [33l. |34|). and second, any approximations in the calculations of the integer- valued replica partition function 
are quite risky for the validity of the further analytic continuation to non-integer values of N. 

For the problem under consideration the point is that, in fact, neglected exponentially small contributions at integers 
N > 1 appear to be quite essential in the region N — ¥ 0, which defines the properties of the free energy distribution 
function P*( f) = hniL-^oo Pl{J)- In other words, the problem is that the two limits, L — > 00 and N — > 0, do not 
commute (47I.|48|. 
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Nevertheless, in terms of this approximation (assuming the universal scaling L 1 / 3 of the free energy fluctuations) 
one can derive the left tail of the distribution function P*(/) which is given by the asymptotics of the Airy function, 

P*{f ~~ * — oo) ~ cxp(— 1|/| 3 2 ) (4t|. For the first time the form of the right tail of this distribution function, 
P*{f ~ ► +00) ~ exp[— (const) f 3 ] , has been derived in terms of the optimal fluctuation approach [50l - l52l ]. where it has 
been demonstrated that both asymptotics (left and right) of the function P*(f) are consistent with the Tracy- Widom 
distribution |l| which was known to describe the statistical properties of many other systems [2|-|7| . 

For the first time TW distribution for the directed polymers with (5-correlated random potentials was derived in 
terms of the distribution of the solutions of the KPZ equation [ij [l(J, which, in particular, describes the domain walls 
growth, and which is equivalent to the present system. Almost simultaneously the exact solution of the one-dimensional 
directed polymer problem has been found in terms of of the replica method, which involved the summation over the 
whole spectrum of excited states in the corresponding JV-particle quantum boson system [liTfl4j |. These calculations 
resulted in the derivation of the entire free energy distribution, which was proved to coincide with TW distribution 
function. 

The Tracy- Widom function, Eq. p.2ip . was originally derived for rather specific mathematical problem, namely for 
the probability distribution of the largest eigenvalue of a JV x JV random hermitian matrix in the limit JV — > 00 [l[ . 
It is amazing but now there are exists a long list of statistical systems (which at first sight have just nothing to do 
with the original random matrix problem) whose macroscopic properties are described by the same universal TW 
distribution function. In other words, all the above observations indicate there exist a kind of "superuniversality" for 
the entire class of various random systems. 

In this review we have described the exact solution of the problem which remained unsolved during last almost 
thirty years. It should be stressed that this solution has been obtained in terms of the replica method. This is 
very rare case when the solution of a non-trivial problem has been found without using heuristic "replica magic" 
operations, quite typical for this method, which usually forced to think that the "replica method" and the "exact 
solution" are two things absolutely incompatible. Hopefully finding exact solution is not always means the end of 
the story: methodology and created mathematical technique could be used for solving numerous other problems still 
waiting for their solutions... 



Appendix A 
Quantum bosons with repulsive interactions 



1. Eigenfunctions 

The eigenstates equation for JV-particle system of one-dimensional quantum bosons with (5-interactions is 

N N 

-J2d 2 x *(x) + ^X;<y(s«-s&)*(x) = -JJ£*(x) (A.l) 

(where k = j3 3 u). Due to the symmetry of the wave function with respect to permutations of its arguments it is 
sufficient to consider it in the sector 

x\ < X2 < ■■• < %n (A.2) 
as well as at its boundary. Inside this sector the wave function ^(x) satisfy the equation 

JV 

2 



1 N 

-£e£*(x) = -025*(x) (A.3) 



which describes JV free particles, and its generic solution is the linear combination of JV plane waves characterized by 
JV momenta {qi, q2, gjv} = q- Integrating Eq. (|A.l|) over the variable (xi+\ — Xi) in a small interval around zero, 
\x i+1 — Xi\ < e —> 0, and assuming that the other variables {xj} (with j 7^ i, i + 1) belong to the sector, Eq. (|A.2[) . one 
easily finds that the wave function Vf'(x) must satisfy the following boundary conditions: 



(d Xi+1 -9 S4 + k)*(x) 



= (A.4) 
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Functions satisfying both Eq. (|A.3[) and the boundary conditions Eq. (|A.4[) can be written in the form 



* ft .„< w (x 1 ,...,xjr) = tfW(x) = C (Y[[d Xa -d Xb +K\) det exp(i 



N 



a<b 



iq C Xdj 



(c,d) I V 



(A.5) 



where C is the normalization constant to be defined later. First of all, it is evident that being the linear combination 
of the plane waves, the above wave function satisfy Eq. ([A.3p . To demonstrate which way this function satisfy the 
boundary conditions, Eq. (|A.4[) . let us check it for the case i = 1. According to Eq. (|A.5j) . the wave function \l/q (x) 
can be represented in the form 



*W(x) = -(a BS -a B1 -«)*W(x) 



where 



#W(x) = C J] - + k] [d X2 - d Xa +«] ) ( J] - a xt + k] det expfe 



a— 3 



JV 



3<a<6 



(c,d)=l,...,iV 



(A.6) 



(A.7) 



One can easily see that this function is antisymmetric with respect to the permutation of x\ and Xi- Substituting 
Ea.(|A.6Jl into Eq. (|A.4l) (with i = 1) we get 



(9 X2 - 9 Xl ) 2 - k 2 



*W(x) 



= 



(A., 



Given the antisymmetry of the l.h.s expression with respect to the permutation of x\ and xi the above condition is 
indeed satisfied at boundary x\ = X2- 

Since the eigenfunction vtq^x) satisfying Eq. (|A.ip must be symmetric with respect to permutations of its argu- 
ments, the function, Eq. (|A.5p . can be easily continued beyond the sector, Eq. (|A.2p . to the entire space of variables 
{x 1 ,x 2 , -,x N } € R N , 



N 



*W(x) = c(]l[-i(d Xa -d Xb )+iKsgn(x a -x b ) 



det 



exp(ig c x d ) 



(c,d) = l,...,N 



(A.9) 



where, by definition, the differential operators d Xa act only on the exponential terms and not on the sgn(x) functions, 
and for further convenience we have redefined i Ar ( Ar_1 )/ 2 C' — > C. Explicitly the determinant in the above equation is 



det 



exp(iq c x d ) 



(c,d)=l,...,JV 



N 

exp[*£, 



(A.10) 



where the summation goes over the permutations P of N momenta {91,92, ■■■tQn} over N particles {xi,X2, ...,a;jv}, 
and [P] denotes the parity of the permutation. In this way the eigenfunction, Eq. (|A.9p . can be represented as follows 



, N v N 

*W(x) = CJ2(-V IP] (l[[-i{d Xa ~d Xb ) + ik sgn(x a -x b )\) exp[z^g Pa x a 



a<b 



Taking the derivatives, we obtain 



/ N v N 

= C^](-1) 1P1 ( ]j[qp a ~q Pb +insgn(x a - x b ) J exp ^ g Pa x a 

P \<b ' a=l 



(A.ll) 



(A.12) 



It is evident from these representations that the eigenfunctions \& q (x) are antisymmetric with respect to permuta- 
tions of the momenta q\,...,qN- 

Finally, substituting the expression for the eigenfunctions, Eq. (|A.5|) (which is valid in the sector, Eq. (|A.2jl ). into 
Eq. (|A.3|) for the energy spectrum we find 



E = 



1 N 



2/3 



(A.13) 
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2. Orthonormality 

Now one can easily prove that the above eigenfunctions with different momenta are orthogonal to each other. Let 

q' 



us consider two wave functions (x) and (x) where it is assumed that 



qi < <72 < ■■■ < QN 
q[ < q' 2 < ... < q' N 



(A.14) 



Using the representation, Eq. (|A.lip , for the overlap of these two function we get 

/+oo 
^x*^ (x)*W(x) 
-OC 



\C\ 2 J2 (-1) [P1+[P ' ] [ + °° d N J (f[ [i(d Xa - d Xb ) ~ i KS gn(x a - x b )]) expH]T q' p , x a ] } 
p,p> J -°° \< b ' 0=1 ' 

N 



H[ _i (^x a -9x b ) + iKSgn(x a - x b )] J exp[i^2q Pa x a ] > 

K a<b ' a=l ' 

Integrating by parts we obtain 

r +oc N 



(A.15) 



(A.16) 



a=l 



A' 



Il[ _ *(^a _ &xb) - *«sgn(x a - x b )] [-i(d Xa ~ d Xb ) + i/csgn(x a - x b )] J exp \i ^ q Pa x a ] 



a<b 



N 



r +oo N , N \ N 

*W*(x)*W(x) = |C| 2 J2(-1) [P]+[P ' ] / ^xexp[-^g;^ a ] nh( 9 - - 9 -) 2 + «'] exp[»X;*.*.] 

P,P' J -°° a=l \i<b ' o=l 



Taking the derivatives and performing the integrations we find 



(A.17) 



/ N \ r + oo N 

*W>)*W(x) = |C| 2 ^(-l)^'Un[(^-fe) 2 + « 2 ]) / /xe X p[^( gPo - g ;,) Ia ] 

P.P' \<b ' J -°° 0=1 



P.P' 



(A.18) 



Taking into account the constraint, Eq. (|A.14[) . one can easily note that the only the terms which survive in the above 
summation over the permutations are P — P', all contributing equal value. Thus, we finally get 



N 



*W'(x)*W(x) = |C| 2 Nl [( 9b - q b f + K 2 ] H(2«)5(q a q' a 



a<b 



r N 



L a=l 



This relation defines the normalization constant 

i2 



C(q) 



^n.< b [(9a- 9 .) 2 + K 2 ] 



(A.19) 



(A.20) 



The proof of completeness of this set is given in Ref. [4l|. It should be noted that the above wave functions present 
the orthonormal set of eigenfunctions of the problem, Eq. (IA.l[) . for any sign of the interactions k, e.i. both for the 
repulsive, k < 0, and for the attractive, k > 0, cases. However, only in the case of repulsion this set is complete, while 
in the case of attractive interactions, k > 0, in addition to the solutions, Eq. (|A.lll) . which describe the continuous 
free particles spectrum, one finds the whole family of discrete bound eigenstates. Detailed description of these states 
is given in Appendix B. 
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Appendix B 
Quantum bosons with attractive interactions 



1. Ground state 



The simplest example of the bound eigenstate defined by eq. flA.ll) (with n > 0) is the one in which all N particles 
are bound into a single " cluster" : 



*W(x) = C exp 



N N 

X a - ^2 \ X a 



Xb\ 



a, 6=1 



(B.l) 



where C is the normalization constant (to be defined below) and q is the continuous momentum of free center of mass 
motion. Substituting this function in Eq. (|A.ip , one can easily check that this is indeed the eigenfunction with the 
energy spectrum given by the relation 



JY 



JY 



E = ~2~B^ iq - ^nJ2 sgn(x a - x b ) 

P a=l <- 



6=1 



(B.2) 



where it is assumed (by definition) that sgn(0) = 0. Since the result of the above summations does not depend on the 
mutual particles positions, for simplicity we can order them according to Eq. (|A.2[ ). Then, using well known relations 



N 



sgn(x a - x b ) = -(JV + l-2o) 
1 



6=1 



A' 



J^a = ^N(N + 1) 

a=l 
N 

E« 2 



o=l 



N 

(N+1)(2N+1) 

6 



for the energy spectrum, Eq. (|B.2p . we get 



N 



The normalization constant C is defined by the orthonormality condition 

/+oo 
dxi...dx N *y (x)*«(x) = (27r)%-</) 
-oo 

Substituting here Eq. (|B.l[) we get 

,.+00 r N N 

*J ) *(x)f( 1) (x) = \C\ 2 / dx x ...dx N exp i{ q - q ')Y J X a --K J2 K 

J — oo L 1 * „ J._i 



(B.3) 
(B.4) 
(B.5) 

(B.6) 
(B.7) 



Xb\ 



a, 6=1 



/•+oo /* + oo 


r+oo 




/ dx\ I dx2--- 


. dxN exp 


i(q 


J —oo J X\ 







N 



N 



l(q - q')J2x a + nJ2(N + 1 - 2a)x a 



(B. 



where for the ordering, Eq. (|A.2[) . we have used the relation 



1 N 

- \%a-x b \ = -^(JV +l-2a)a; a 



(B.9) 



a, 6=1 



a=l 
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Integrating first over xn, then over xn-i, and proceeding until x\, we find 



^(x^x) = |C| 2 AT! J] 



^ r[(iV — r)n — i(q — q')] 
N-l 



dxi exp [iN(q — q')x\\ 



1 



\^{U r{N _ r)K 

|C|2 lV!^ (2-)^-^) 



(2ir)5{N(q-q')) 



According to Eq. (|B.7p this defines the normalization constant 

C = 



i N Nl 



kN 



(B.IO) 



(B.H) 



Note that the eigenstate described by the considered wave function, Eq. (|B.lj) . exists only in the case of attraction, 
K > 0, otherwise this function is divergent at infinity and consequently it is not normalizable. 



It should be noted that the wave function, Eq. (|B.I[) , can also be derived from the general eigenfunctions structure, 
Eq. (|A.I2j) . by introducing (discrete) imaginary parts for the momenta q a . We assume again that the position of 
particles are ordered according to Eq. (|A.2j ). and define the particles momenta according to the rule 



q a = q-- K (N + l-2a) 

Substituting this into Eq. (|A.I2[) we get 

< < < X N ) OC )[ p l ( [] 

P ^a<b 

r N N 

iqJ2 X a + ~ J2( N + 1 - 2Pa)Xa 



(B.I2) 



N r 



q - -k(N + 1 - 2P a ) 



q - - K (N + 1 - 2P b ] 



x exp 



a=l 



/ N s r N N 

E(- 1 ) [P1 II - p « + !] ) ex p x * + 1 J2( N + 1 - 2P *>* 



a<b 



a=l 



(B.I3) 



Here one can easily note that due to the presence of the product Il^<6 [Pb ~ Pa + 1] in the summation over permutations 
only the trivial one, P a = a, gives non-zero contribution (if we permute any two numbers in the sequence 1,2, N 
then we can always find two numbers a < b, such that P^ = P a — 1). Thus 



tfj^fci <x 2 < ... < x N ) 



oc exp 



N ^ N 



a=l 



a=l 



(B.14) 



Taking into account the relation, Eq. (|B.9[ ). we recover the function, Eq. (|B.l[) . which is symmetric with respect to its 
N arguments and therefore can be extended beyond the sector, Eq. (|A.2l) . for arbitrary particles positions. Finally, 
substituting the momenta, Eq. (|B.12p . into the general expression for the energy spectrum, Eq. (|A.13j) . we get 



E = YB^ q -r [N+l - 2a) ^ 

^ a=l 



N 



(B.15) 



Performing here simple summations (using Eqs. (|B.4p . (|B.5|) ) one recovers Eq. (|B.6|) . 



2. Eigenfunctions 
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A generic eigenfunction of attractive bosons is characterized by N momenta parameters {q a } (a = 1, 2, ...N) which 
may have imaginary parts. It is convenient to group these parameters into M (1 < M < N) "vector" momenta, 

q? = q a - ~ « (n a + 1 - 2r) (B.16) 

where q a (a = 1,2, ...,M) are the continuous (real) parameters, and the discrete imaginary components of each 
"vector" are labeled by an index r = 1, 2, n a . With the given total number of particles equal to N, the integers n a 
have to satisfy the constraint 

M 

n a = N (B.17) 

In other words, a generic eigenstate is characterized by the discrete number M of complex "vector" momenta, by the 
set of M integer parameters {n\ 1 ri2, tlm } = n (which are the numbers of imaginary components of each "vector") 
and by the set of M real continuous momenta {qi, q?, <zm} = q- 

The general expression for the eigenfunctions is given in Eqs. (|A.9P - (|A.12p . To understand the structure of the 
determinant of the N x N matrix exp(iq a xi,) , which dchncs these wave functions, the N momenta q ai eq. (|B.16p . can 
be ordered as follows: 



{qa} = {q?} = {q{, ql, - , <&; ql ql - , ql 2 ] - ; qf , q% - . <7^t} 



(B.18) 



By definition, 



det 



exp(ig a x c ) 



(c,d)=l,...,N 



N 

^(-l) [p] exp\i^2q Pa x a 



(B.19) 



where the summation goes over the permutations of N momenta {q a }, Eq. (|B.18p . over N particles {x\, X2, xn}, and 
[P] denotes the parity of the permutation. For a given permutation P a particle number a is attributed a momentum 

component The particles getting the momenta with the same a (having the same real part q a ) will be called 

belonging to a cluster For a given permutation P the particles belonging to the same cluster are numbered by 
the "internal" index r = 1, ...,n a . Thus, according to Eq. (|A.lip . 



N 



Yi [~*(^a - 9 Xb ) + insgn(x a - x b ) 

a<b 



exp 



N 



a(c) 
(c) ' 



(B.20) 



where Cq jn is the normalization constant to be defined later. Substituting here Eq. (|B.16l) and taking derivatives we 
get 



N r 



= C^^-^U (?«(a)-y [na(a) + l-2r(a)])-(g a(6 )-y[n a(6) +l-2r(6)])+i«Bgn(x -a:6) 



a<b 



x exp 



i X! &*(c)Z c + ^ X! (""(c) + 1 _ 2r ( c )) : 



(B.21) 



The pre-exponential product in the above equation contains two types of term: the pairs of points (a, b) which belong 
to different clusters (a(a) =^ &(b)), and pairs of points which belong to the same cluster (a (a) = a(b)). In the last 
case, the product n a over the pairs of points which belong to a cluster fi Q reduces to 



Ila oc Yl [ r ( b ) ~ r ( a ) ~ s S^( x a ~ Xb)}] 



(B.22) 



a<beO Q 



Similarly to the ground state wave function Eq. (|B.13[) (|B.14I) . one can easily note that due to the presence of this 
product in the summations over n Q ! "internal" (inside the cluster Q a ) permutations r(a) only one permutation gives 
non-zero contribution. To prove this statement, we note that the wave function t&qjn (x) is symmetric with respect 
to permutations of its N arguments {x a }; it is then sufficient to consider the case where the positions of the particles 
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are ordered, x\ < x 2 < ■ ■ ■ < xpj. In particular, the particles {x ak } (Ar = 1,2,..., n a ) belonging to the same cluster 

In this case 



£l a are also ordered x ai < x a2 < ■ ■ ■ < x 



h„ ex JJ[ r (0-r(fc) + l] 



(B.23) 



k<l 



Now it is evident that the above product is non-zero only for the trivial permutation, r(k) = k (since if we permute 
any two numbers in the sequence 1, 2, ...,n a , we can always find two numbers k < I, such that r(l) = r(k) — 1). In 
this case 



n Q oc JJ[J-fr+i] 



(B.24) 



k<l 



Including the values of all these "internal" products, Eq. (IB.24p , into the redefined normalization constant Cq^ , for 
the wave function, Eq. (|B.21j) (with x\ < X2 < • • ■ < xjy), we obtain 

Nr. 

*^(x) - ^)^'(-l) [P1 LI («»(«) - y^)) " - +«(Ka)-r(6) - 1)) 



x exp 



a<6 
a(a)^a(fc) 

AT 



AT N 
i E 1ot{c)Xc + ~ E ( n a(c) + 1 ~ 2r (c) 



(B.25) 



where the product goes only over the pairs of particles belonging to different clusters, and the symbol J2p' means 
that the summation goes only over the permutations P in which the "internal" indexes r(a) are ordered inside each 
cluster. 

Now taking into account the symmetry of the wave function "f^n (x) with respect to the permutations of its 
arguments the expression in Eq. (|B.25p can be easily continued beyond the the sector X\ < x 2 < ... < xjy for the entire 
coordinate space Rn- Using the relations 



(n a + 1 - 2r(a))x a = ^(n Q + 1 - 2k)x ak = -- ^ \x ak 



(B.26) 



a£fl a 



fc=l 



k,l=l 



(where x ai < x a2 < ... < x ana ), for the wave function "fq^n (x) with arbitrary particles positions we get the following 
sufficiently compact representation (cf. Eq. (|B.20p '): 



JV 



cw^s , E(- i ) |pi n 



exp 



M n a M n a 

E *<* E Xc - 7 E E 



-!(^ ~d Xh ) + insga{x a - x b ) 

P a<b 

(B.27) 

Note that although the positions of particles belonging to the same cluster are ordered, the mutual positions of 
particles belonging to different clusters could be arbitrary, so that geometrically the clusters are free to "penetrate" 
each other. In other words, the name "cluster" does not assume geometrically compact particles positions. 
Finally, substituting Eq. (IB.16[) - (IB.17[) into Eq. (|A.13[) . for the energy spectrum one easily obtains: 



E 



M 



2/3 ^ £^™> 2/3 

a — 1 r=l a— 1 a — 1 



(B.28) 



3. Ortho normality 

We define the overlap of two wave functions characterized by two sets of parameters, (M, n, q) and (M', n', q') as 

f+co 

Qiy !) M)^ / d"**£J(*)*<# (x) (B.29) 

J — CO 
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Substituting here Eq. (IB.27[) we get 

, , Z' + OO 

O W) = C , C"EE(- 1 ) im|P1 / d 

r> r>/ J — CO 



(B.30) 



p P' 



N 

n 

a<b 
a (a)^a' (b) 



*(<9x a ~d Xb ) - iKSgn(x a - x b ) 



N 

n 

a<b 
a(a)^a(b) 



-i(d Xa - d Xb ) +i«;sgn(x a - x b ) 



exp 



exp 



M' "a M' «4 

-* E ^ Y Xc - 1 E Y \ Xc ' 



a=l c££2' 



A I 



n a M n a 

1 Y q<* Y Xc - 1 Y Y \ x c-xc 



a=l c,c'GO Q 



where {f2 a } and {^' Q } denote the clusters of the permutations P and P' correspondingly. Integrating by parts we 
obtain 



qLT W) = c^cffJ E E M) [p]+[p,] / < 

r> r>r '' — CO 



d N 



x exp 



p p> 



Y <t« Y Xc - 7 Y Y \ Xc 



N 

n 

a<b 
a' (a)^a' (b) 



-i(d Xa - d Xb ) - iKSgn(a; Q - x b ) 



a=l ceiV a 
N 



a—1 c,c'6£i' 



n 

a(a)^a(b) 



-i(d Xa - d Xb ) + iKSgn(x a - x b ) 



x exp 



M 



E qa E Xc _ 7 E E l x c - a; c 



a=l cef2 Q 



a— 1 c,c'Efi a 



(B.31) 



First, let us consider the case when the integer parameters of the two functions coincide, M = M', n = n', and 
for the moment let us suppose that all these integer parameters {n a } are different, 1 < m < < ... < Um- Then, 
in the summations over the permutations in Eq. (IB.31|) . we find two types of terms: 

(A) the "diagonal" ones in which the two permutations coincide, P = P' ; 

(B) the "off-diagonal" ones in which the two permutations are different, P ^ P' . 
The contribution of the " diagonal" ones reeds 



Qi M n M)(A) (q, q') = c^c^r e / «P Y £ Y »- - \ Y Y 

p , '~ 00 <- a=l ceQ a a=l cc'ef 



Xr Xr 



N 

n 

a<b 



(d Xa - d Xb ) 2 + k 2 



exp 



M n a M n a 

Y^Y^-jY Y 



a—1 c,c'££] Q 



(B.32) 



It is evident that all permutations a{a) in the above equation give the same contribution and therefore it is sufficient 
to consider only the contribution of the "trivial" permutation which is represented by Eq. (|B.18j) . The cluster ordering 
given by this permutation we denote by ao(a). For this particular configuration of clusters we can redefine the 
particles numbering, so that instead of a "plane" index a — 1, 2, N the particles would be counted by two indexes 
(a, r): {x a } — > {x" } (a — 1, M) (r = I, n a ) indicating to which cluster a a given particle belongs and what is 
its "internal" cluster number r. Due to the symmetry of the integrated expression in Eq. (|B.32j) with respect to the 
permutations of the particles inside the clusters, we can introduce the "internal" particles ordering for every cluster: 
x? < x% < •■• < x" .In this way, using the relation, Eq. ([B.26[) . we get 



ni!n 2 !...7j M ! 



M 

n 



dxl 



4- co 



+ CO 



dxl 



x exp 



M n a 



a—1 r—1 
M n a np . 



M n a 



a=l r=l 



E^E< + ?EE(^ + 1 - 2r )< 

i e q a Y < + 7 E E( n « + 1 - 2r > 



2 K 2 



nnn -(^-9 < y+ 

a<j3 r—1 r' — 1 L 



exp 



a—1 r—1 



M n a 



(B.33) 
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where the factor N\/ni\...n^ is the total number of permutations of M clusters over N particles. Taking the 
derivatives and reorganizing the terms we obtain 



« m)<a W) = «»c m n n n 



M n a np 



«</3 r—1 r' — l 



{<l a ~ TT™a) - (<7/3 - ir n p) +in{r -r' - 1) 



f /* + 0O /* + 0O /* + 0O I 

„,_i J —00 J x? J x°: 1 



Simple integrations over yields (cf. Eqs. (|B.8|l - (|B.10p ) 

n (M,M)< A » 



/• M n a 

W) = i^rwnnn 

V</3 r=l r' = l 



(Q0 - y"-/?) + iK ( r ~ r ' ~ !) 



M 

n 



(n a !) 2 K n ° 



(27T)5(q a - q' a ) 



(B.35) 



Now let us prove that the "off-diagonal" terms of Eq. (|B.31l) . in which the permutations P and P 1 are different, 
give no contribution. Here we can also chose one of the permutations, say the permutation P, to be the "trivial" one 
represented by Eq. (|B.18|) with the cluster ordering denoted by 0.0(0). Given the symmetry of the wave functions it 
will be sufficient to consider the contribution of the sector x\ < X2 < ■■■ < xn- According to Eq. (|B.3ip . we get 



„ ,- M n a M n a 

Q(f n ^ <B, (q,q') OC ^(-1)^ / d?* -p -iJ2<l' a E ^-lE E I*. 

P' Jx 1 <...<x N L Q=1 Qen , a=la.beO'„ 



Xb\ 



N 

n 

a<b 
a (a)^a (b) 



-i^xa - dxb) + «KSgn(ir a - x b ) 



N 

n 

a<b 
a' (a)^a/ (b) 



-i(dx a -9 Xi )- in sgn(x a - x b ) 



x exp 



M n a 



M n a 



' E q ° E x °- + ? E E ( n ° + 1 ~~ 2r ( a )) X a 



a—1 a£0° 



a — 1 a^n° 



(B.36) 



Here the symbols } denote the clusters of the trivial permutation ao(a). Since P' 7^ P, some of the clusters D,' a 
must be different from f2° . As an illustration, let us consider a particular case of -/V = 10, with three clusters rii = 5 
(denoted by the symbol "O") 1 n 2 = 2 (denoted by the symbol " x") and = 3 (denoted by the symbol "A"): 



particle number a 


1 


2 


3 


4 


5 


G 


7 


8 


9 


10 


permutation ao(a) 

















X 


X 


A 


A 


A 


permutation a' (a) 











A 





X 


X 





A 


A 



Here in the permutation a' (a) the particle a = 4 belong to the cluster a — 3 (and not to the cluster a — 1 as in 
the permutation Qo(a)), and the particle a = 8 belong to the cluster a = 1 (and not to the cluster a = 3 as in the 
permutation ato(a)). Now let us look carefully at the structure of the products in Eq. (|B.36[) . Unlike the first product, 
which contains no "internal" products among particles belonging to the cluster f2J, the second product does. Besides, 
the signs of the differential operators (d Xa — d Xb ) in the second product is opposite to the "normal" ones in the first 
product (cf. Eqs. (|B.22j) - (|B.24|) '). It is these two factors (the presence of the "internal" products and the "wrong" 
signs of the differential operators) which makes the "off-diagonal" contributions, Eq. (IB.36[) . to be zero. Indeed, in the 
above example, the second product contains the term 



114,5 = 



i{d x 



exp 



3 n a 

z E^ E 

a—1 a€Q% 



\ E E ("<* + 1 - 2r (a)K 

a—1 aGO° 



(B.37) 
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(we remind that the particles in the clusters f2° are ordered, and in particular 24 < x$) 
get 



Taking the derivatives, we 



114,5 



- (kx + I (ni + 1 - 2r(4)) - % - | (rn + 1 - 2r(5)) j + k 



cx [r(4) - r(5) + l] 







(B.38) 



since in the first cluster r(a) — a. 

One can easily understand that the above example reflect the general situation. Since all the cluster sizes n a are 
supposed to be different, whatever the permutation a'{a) is, we can always find a cluster f2° such that some of its 
particles belong to the same cluster number a in the permutation a' (a) while the others do not. Then one has to 
consider the contribution of the product of two neighboring number points 



n' 



k,k+l 



-i{d Xk -d Xk+1 ) + in 



exp 



M 



M 
K \ ~\ 



J~] (n a + 1 - 2r(a))x a 



a=l o6!l° 



(B.39) 



where in the permutation a' (a) the particle k belong to the cluster number a and the particle {k ■ 
other cluster. Taking the derivatives one gets 



ni 







1) belong to some 
(B.40) 



ife, fe +i « [Kfc) - r(k + 1) + 1] 
as r(a) is the "internal" particle number in the cluster fi°, where r(k + 1) = r(k) + 1 (cf. Eqs. (|B.22l) - (|B.24p ). 

Thus, the only non-zero contribution to the overlap, Eq. (|B.29p . of two wave function \&^2( x ) and "J^n (x) (having 
the same number of clusters M and characterized by the same set of the integer parameters 1 < ri\ < n% < ... < um) 
comes from the "diagonal" terms, Eq. (|B.35j) : 



M 



)(iW,M) 
' n,n 



(q,q') 



\ C {M)\ 
I q.n I 



M 



MR 



Tify Hi 



]\2 K n a 



M n a rip 

nnn 

a</3 r=l r' = l 



Via - -7rn a - (qp 



-np) + in (r — r 1 — 1) 



n 



(2n)S(q a - q' a ) 



(B.41) 



The situation when there are clusters which have the same numbers of particles n a is somewhat more complicated. 



Let us consider the overlap between two wave function ^^f^i^) and ^^■'(x) (which, as before have the same M 



and n) such that in the set of M integers rii,n 2 , ...,nM there are two n a 's which are equal, say n ai 



(where 



Oil 7^ a 2)- In the eigenstate (q',n) these two clusters have the center of mass momenta q' ai and q' a2 , and in the the 
eigenstate (q, n) they have the momenta q ai and q a2 correspondingly. According to the above discussion, the non-zero 
contributions in the summation over the cluster permutations a{a) and a' (a) in Ea. (|B.31[) appears only if the clusters 
{f2 Q } of the permutation a{a) totally coincide with the clusters {£l' a } of the permutation a'(a). In the case when 
all n a are different this is possible only if the permutation a (a) coincides with the permutation a' (a). In contrast to 
that, in the case when we have n ai = n a2 , there are two non-zero options. The first one, as before, is given by the 
"diagonal" terms with a(a) = a'(a) (so that the clusters {fl a } and {Q' a } are just the same), and this contribution 
is proportional to 5(q ai — q' ai )5{q a2 — q' a2 )- The second ("off-diagonal") contribution is given by such permutation 
a' (a) in which the cluster Q' ai (of the permutation a' (a)) coincide with the cluster fi a2 (of the permutation a(a)) and 
the cluster £l' a2 (of the permutation a' (a)) coincide with the cluster J7 Ql (of the permutation a(a)) while the rest of 
the clusters of these two permutations are the same, fl' a = Q a (a 7^ cti,ct2)- Correspondingly, this last contribution 
is proportional to 5{q ai — q' a2 )S(q a2 — q' ai ) (— l) n °i. In fact this situation with two equivalent contributions is the 

consequence of the symmetry of the wave function ^ / l^n( x ) : the permutation of two momenta q ai and q a2 belonging 



to the clusters which have the same numbers of particles, n ai 



produces the factor (— l) n °i. This is evident 



from the general expression for the wave function, eq. (IA.9l) . where the permutation of any two momenta q ai and q a2 
belonging to the clusters which have the same numbers of particles corresponds to the permutation of n columns of 
the matrix exp(iq Q Xb) . Therefore considering the clusters with equal numbers of particles as equivalent and restricting 
analysis to the sectors q ai < q a2 ; q' ai < q' a2 we find that the second contribution, 8(q ai —q' a2 ) 5(q a2 —q' ai ) is identically 
equal to zero, thus returning to the above result Eq. (|B.41[) . 

A generic eigenstate (q, n) with M clusters could be specified in terms of the following set of parameters: 



(q, n ) = {(<?!> mi),...,(« Sl , mi); (q Sl+1 ,m 2 ), ■-, (q Sl +s 2 ,m 2 ); .... ; (q Sl+ ... +Sk _ 1+1 ,m k ), (q Sl 



+ ... + s fc 



TOfc)} (B.42) 
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where s\ + S2 + ■■■ + Sk = M and k integers {mj (1 < k < M) are all supposed to be different: 

1 < mi < 7712 < ■•■ < fnk (B.43) 
Here the integer parameter k denotes the number of different cluster types. For a given k 



M 



(B.44) 



a=l 



Due to the symmetry with respect to the momenta permutations inside the subsets of equal n's it is sufficient to 
consider the wave functions in the sectors 



qi < q 2 < ... < q Sl ; 

9si+i < q Sl +2 < •■• < q Sl +s 2 ; 



(B.45) 



&! + ... +«*-!+! < <?Sl + ...+Sfc_l+2 < •■• < q Sl +...+S k - 1 +S 



In this representation we again recover the above result Eq. (IB.4ip 

Finally, let us consider the overlap of two eigenstates described by two different sets of integer parameters, n' 7^ n. 
In fact this situation is quite simple because if the clusters of the two states arc different from each other, it means 
that in the summation over the pairs of permutations P and P' in Eq. (IB.3ip there exist no two permutations for 
which these two sets of clusters {f2 Q } and {£l' a } would coincide. Which, according to the above analysis, means that 
this expression is equal to zero. Note that the condition M 1 ^ M automatically implies that n' 7^ n. 

Thus we have proved that 



Q 



(M,M') , t \ 



q.q') = K^f S(M,M') J] S(n a ,n' a ) Y[(27r)S(q a q' c 



M 



M 



dvi n 



(»„!) 2 k"» 



M n a rip 

nnn|(«- 

a<8 r=l r' = l 



yn a ) - (qp - —np) + in [r - r' - 1) (B.46) 



where the integer parameters {n a } and {n' a } are assumed to have the generic structure represented in Eqs. (IB.42|) - 
(|B.44jl . and the momenta {q a } and {q' a } of the clusters with equal numbers of particles are restricted in the sectors, 
Eq. (|B.45p . According to Eq. (|B.46l) . the orthonormality condition defines the normalization constant 



|C(-) (q ,n)| 2 = ^ 



r M 

n 

a=l 



{n a \) 2 K n ° 



M n a Up 

nnnr — — — — ^— 

a</3 r=l r'=i (q a - fn a ) - (q p - fnp) +in(r-r' - 1) 



(B.47) 



4. Propagator 

The time dependent solution ^(x, t) of the imaginary-time Schrodinger equation 

N „ N 



with the initial condition 



*(x;0)=C =1 5(x a ) 



can be represented in terms of the linear combination of the eigenfunctions ^^/(x), Eq. (IB. 27ft : 



*(x,t) = J2 E' /'^q *^( x )*^"(0) ex P [-£ M (q,n)i] 



(B.48) 



(B.49) 



(B.50) 
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where the energy spectrum Em(q, n) is given by Eq. (|B.28[) . The summations over n a are performed here in terms of 
the parameters {s l ,m i \, Eqs. (|B.42|> - (jB.44j) : 



E=E E E *(£-<. 

n fc=l si...sj.=l l<mi...<m fc M=l ' M=l ' 



(B.51) 



where 8(n,l) is the Kronecker symbol, and for simplicity (due to the presence of these Kronecker symbols) the 
summations over rrii and are extended to infinity. The symbol f T>q in Eq. (IB.50|) denotes the integration over M 
momenta q a in the sectors, Eq. (|B.45l) . 

The replica partition function Z(N, L) of the original directed polymer problem is obtained via a particular choice 
of the final-point coordinates, 



i r' 

Z(N,L) = #(0;L) = E E / Vc ± 1*^(0)1 exp[-E M (q, n) L] 

71/T 1 n " 

According to Eq. (lB~27l) . for M > 2, 



(B.52) 



M=l n 



AT 

= ^E'(-i) [p] n 



a<b 
*(a)^a(i>) 



9a(a) - ylf. 



a(a) + 1 - 2r ( a )] ) ~ («a(&) ~ y [«a(&) + 1 - 2r ( & )] ) 



TV! 



qn m!n 2 !...n M ! 



M n a r 

nnn 



9a - TT n a 1-19/3- V n /3 + ZK 



a</3 r=l r'=l 

Substituting here the value of the normalization constant, Eq. (|B.47jl . we get 



(B.53) 



|«(o)| 2 = 



NIk 



v g Q - - (q^ - fnp) + in (r - r' 



(riaii «n a ) a <p n"=x n"li (9a - f n a ) - ^ - fn p ) + 1* (r - V - i; 



(B.54) 



This expression can be essentially simplified. Shifting the product over r' in the denominator by 1 we obtain 



|«(0) 



N\k 



N 



nf=i Kn * ) a<p Ur=i («<* - fn a ) - (qp - fnp) +in(r-n p - 1) 



M 

n 



n"=i (9« - f n„) - - f n^) + we (r - 1) 



(B.55) 



Redefining the product parameter r in the denominator, r — > n a + 1 — r, and changing the obtained expression (under 
the modulus square) by its complex conjugate we get 



|^(o)| 2 = 



Wk n " llr=i (9a - f n a ) - (q p - fnp) + ik (r - 1) 



n 



(riali Kn «) <*<fi n"=i (fc - f n a ) - (g/a - f np) + us r 



(B.56) 



Shifting now the product over r in the numerator by 1 we finally obtain 



1*^(0) 



N\k n 



M 

n 



\q a - qp - t ( n a - np)\ 



(B.57) 



For M = 1, according to Eqs. (lB~T1) and ([HTTT]> . 



|*< 1} (o)| 



(B.58) 



Since the function /(q, n) = ^q^CO)] 2 exp[— Sm(Q; n ) L\ m Eq. (|B.52l) is symmetric with respect to permutations 
of all its M pairs of arguments (q a , n a ) the integrations over M momenta q a can be extended beyond the sector defined 
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in Eq. (|B.45[) for the whole space Rm- As a consequence, there is no need to distinguish equal and different n Q 's any 
more, and instead of Eq. (|B.5f I) , we can sum over M integer parameters n a with the only constrain, Eq. ([B.17p (note 
that this kind of simplifications holds only for the specific "zero final-point" object ^>(0;t), Eq. (|B.52j) . and not for the 
general propagator ^(x; t), Eq. (|B.50p containing N arbitrary coordinates x%, xjy). Thus, instead of Eq. (|B.52p we 
get 



N r M r+oa , oo -, / M x 

WL) = E i n / -t E 5 (E^> A l«(o)| 2 e 



-E M (q,n)L 



(B.59) 



Substituting here Eqs. (|B.28p . (|B.57|) and (|B.58p we get the following sufficiently compact representation for the replica 
partition function: 



Z(N.L) = N\k 



N 



dq 
2-kk.N 



exp 



N r A/ oo „ + oc 

EinE 

M—1 ' L a=l n a = l 



dq a 



2iTKn a 



M 



A/ 



1 9a - 9/3 - iffaa - n/3)| 



*</3 1 3a - 9/3 - l f(n a +np)\ 



A I 



x exp 



K a=l f a=l 



(B.60) 



The first term in the above expression is the contribution of the ground state (M = 1), while the next terms (M > 2) 
are the contributions of the rest of the energy spectrum. After simple algebra the above replica partition function 
can be represented as follows: 



z{n,l) =e 



Z(N,X) 



where fo — ^(3 4 u 2 — ln(/3 3 w), and 



(B.61) 



Z(N.L) 



dq 



2itkN 



N r M OO 

e n e 

M=2 ' l q=1 n Q = l 



x exp 



L 



M 





K 2 L 


' 24/3' 


dq a 


/ M 

<5 




2t M 

— Vn 8 l 
24/3 ^ Q J 

a = l 


} 



2v d 



1 9a ~ 9/3 ~ T ("a ~ 

Q </3 1 9a - 9/3 - f(n a +np)\ 



(B.62) 



Appendix C 
The Airy function integral relations 

The Airy function Ai(x) is the solution of the differential equation 

y"(x) = xy(x) (C.l) 
with the boundary condition y(x — > +oo) =0. At x — > +oo this function goes to zero exponentially fast 

*(*-+«o* v^7i-p("r 3/2 ) (c - 2) 

while at x — > — oo it oscillates and decays much more slowly: 
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The Airy function can also be represented in the integral form: 

Ai(x) = / ^- expf \z 3 - zx] (C.4) 
J c 2m \3 / 

where the integration path in the complex plane starts at a point at infinity with the argument — 7r/2 < Qt—) < —tt/3 
and ends up at a point at infinity with the argument tt/3 < Q{+) < tt/2. Choosing the argument of the staring point 
= — tt/2 + e and that of the ending point 9r + ) = tt/2 — e where the positive parameter e — > is introduced just to 
provide the convergence of the integration, the integration path in eq. (|C.4|) can be chosen to be coinciding with the 
imaginary axes z = iy. 

Just like in the well known Hubbard- Stratonovich transformation the Gaussian function is used to linearize quadratic 
expressions in the exponential, 

cx P (^ 2 ) = | + °°-^ exp(-ix 2 ) exp(Fx) (C.5) 

the Airy function can be used to linearize the cubic exponential terms: 

exp(-F 3 ) = / dx M{x) exp(Fx) (C.6) 

where the quantity F in eq. (|C.6| is assumed to be non negative. One can easily prove this relation using integral 
representation, eq. (|C.4l) . in which the integration path coincides with the imaginary axes, z — iy, and the quantity F 
is first taken to be pure imaginary, F — > %F . In this case the integration over x results in the factor S(F — yj. Further 
trivial integration over y yields the result exp(— iF 3 /3). Performing the analytic continuation of this expression 
F — > —iF one get the relation eq. (|C.6[) . 

In this appendix two other integral relations with the Airy function will be proved. Namely: 

/+oo 

dp Ai(p 2 +uj 1 +lj 2 ) exp[«p(wi -w 2 )] = 2 2/3 7rAi(2 1/3 a;i) Ai(2 1/3 w 2 ) (C.7) 
-co 

and 

h = dy Ai (y + LJ! Ai y + u; 2 = — ^ ^ 1 U 1 ' (C.8) 

JO w l ~~ w 2 



Using the integral representation of the Airy function, eq. (|C.4p . we get 

r+°° [ dz (I \ 
I\ = / dp - — expl —z 3 — p 2 z — ujiz — L02Z + ipuii — ipu>2 ) (C.9) 

J-oo JC 2m V3 ' 

Denoting z = Z\ and ip = z 2 the above integrals can be represented as follows 

h = 2?r I j ^2 exp(iz 3 + z 1 z 2 2 -a; 1 (z 1 -z 2 )-a; 2 (z 1 +z 2 )) (CIO) 



{2mf "V3 

where the integration path C coincides with the imaginary axes. Introducing new integration variables, 

Z1-Z2 = £, 

ZX+Z2 = V (Cll) 

we obtain 

Redefining, £ — > 2 1 / 3 £ and 77 — > 2 1 ^ 3 rj we finally get 

h = ^ J c §~ ex P (^ 3 -2 1 /3 Wl £) J ^ exp (^3„ 2 i/3 W2 ^ = 2 2 /3 7rA i(2 1 /3 Wl )Ai(2 1 /3 a ; 2 ) (C.13) 
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which proves the relation, eq. (|C.7p . 

To prove the relation, eq. (|C.8p . it is sufficient to take into account that the Airy function satisfies the differential 
equation (jC.lj) . Substituting 



Ai(y + uj) = —— Ai"(y + u) 

y + w 



(C.14) 



into eq.flgUl wc find 



h = dy 



1 



Ai"(2/ + o;i) Ai" (y + uj 2 ) 



1 



^1 — W 2 7o 

1 



dy 



dy 



1 



1 



y + Lu 2 y + ^i 



Ai" (y + uji) Ai" (y + cj 2 ) 



Ai" (y + uj 2 ) a • 



y + ^2 



Ai"(i/ + «i) 



Ai"(i/ + Wl ) AJ 



Ai (y + w 2 ) 



(C.15) 



Substituting here Ai"(y + w 2 ) = (y + w 2 ) Ai(y + oj 2 ) for the first term and Ai" (y + lui) — (y + wi) Ai(y + wi) for the 
second term we get 



1 



U)\ - oj 2 

Simple integration by parts eventually yields 



dy Ai" (y + loi) Ai(y + lo 2 ) - / dy Ai" (y + lo 2 ) Ai(y + u>i) 



Wi - cj 2 



- Ai'(wi) Ai(w 2 ) + Ai'(w 2 ) Ai(wi) 



(C.16) 



(C.17) 



which proves the relation, eq. (|C.8p . 



Appendix D 

Fredholm determinant with the Airy kernel and the Tracy- Widom distribution 



In simplified terms the Fredholm determinant det(l — K) can be defined as follows (for the strict mathematical 
definition see e.g. [461]): 



det(l-if) = 1 + J2 



(-iy 



n=l 



dt\dt 2 ...dt n det 



(D.l) 



where the kernel K = K(t,t') is a function of two variables defined in a region a < (t,t') < b. Equivalently the 
Fredholm determinant can also be represented in the exponential form 



det(l - K) = exp 



OO ^ 

— Tr K r ' 



where 



TrK n = 



dtidt 2 ...dt n K{tx,t 2 )K{t 2 , t 3 )...K(t n , t x ) 



(D.2) 



(D.3) 



In this Appendix the original derivation of Tracy and Widom [1] will be repeated in simple terms to demonstrate 
that the function F 2 (s) defined as the Fredholm determinant with the Airy kernel can be expressed in terms of the 
solution of the Panleve II differential equation, namely 



F 2 (s) = det [l - K A ] = exp 



dt(t-s)q 2 (t) 



(D.4) 



where K A is the Airy kernel defined on semi-infinite interval [s, oo): 

K A (t u t 2 ) = Wi)Mfa)-M(h)Ai(t 2 ) 

h — t 2 

and the function q{t) is the solution of the Panleve II differential equation, 

q" =tq + 2q 3 

with the boundary condition, q(t — > +00) ~ Ai(i). 
Let us introduce a new function R(t) such that 

/•OO 

F 2 (s) = exp - / dtR(t) 

or, according to the definition, Eq. (|D.4l) . 

R(s) = A In det(l-# A ) 
Here the logarithm of the determinant can be expressed in terms of the trace: 

In 



det(l-i^) = -^-Trk\ 

71=1 

dti dt 2 ... dt n K A {t u t 2 )K A (t 2 ,t 3 )...K A (t n ,t 1 ) 



71 = 1 



Taking derivative of this expression we gets 



R(s) = - dt(l-K A ) (s,t)K A (t,s) 

J S 

00 poo pc 

= -K A {s,s)-J2 / dt l / 

n J S J S 



dt 9 



n=2 • 



dt„_i K A (s,ti)K A (ti,t 2 ) ...K A {tn-i,s) 



Substituting here the integral representation of the Airy kernel, Eq. (|D.5|) . 

y* 00 

K A (h,t 2 ) = / dzAi(t! + z) Ai(t 2 + z) 
Jo 



after some efforts in simple algebra one gets 



R(s) 



OO f* OO 



dh / dt 2 Ai(h) (1-K A ) (t u t 2 ) Ai(ta) 



(D. 



(D. 



(D. 



(D. 



(D. 



(D.l 



(D.f 



(D.l 



Taking the derivative of this expression and applying some more efforts in slightly more complicated algebra, 
obtain 



ds 



R(s) = -q 2 (s) 



where 



According to Eq. (|D.I3l) . 



dt(l-K A ) (s,t) Ai(t) 



R(s) = / dtq 2 (t) 



(D.l 



(D.l 



(DJ 
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Let us introduce two more functions 

/>oc />oo 

v(s) = / dii / cft 2 Ai(ii)(l-^A) _1 (ti,*2) Ai'(i a ) (D.16) 



p(«) = J dt(l-K A ) \s,t) Ai'(t) (D.17) 

Taking derivatives of the above three functions g(s), u(s) and p(s), Eqs. (|D.14[) . (|D.16|) and (|D. 17|i . after somewhat 
painful algebra one finds the following three relations: 

q = p - Rq (D.18) 
p' = sq - pR - 2qv (D.19) 
v = -pq (D.20) 

Taking derivative of the combination (R 2 — 2v) and using Eqs. (|D.13p and (|D.20|) . we get 

^-(R 2 -2v) = 2q(p - Rq) (D.21) 

On the other hand, multiplying Eq. (|D.18|) by 2q we find 

^-q 2 =2q(p - Rq) (D.22) 

Comparing Eqs. (jD.2ip and (|D.22[) and taking into account that the value of all the above functions at s — > oo is zero, 
we obtain the following relation 

R 2 - 2v = q 2 (D.23) 

Finally, taking the derivative of Eq. (|D.18j) and using Eqs.dDTU, (jDT8|, ((DTT9|> and (fDT23|) we easily find 

q" = 2q 3 + sq (D.24) 

which is the special case of the Panleve II differential equation [37], IH, UlI ■ Thus, substituting Eq. (|D.15|) into Eq. (|D.7[) 
we obtain Eq. (|D.4|) . 

In the limit s oo the function q(s), according to its definition, eq. (|D.14j) . must go to zero, and in this case 
Eq. (|D.24[) turns into the Airy function equation, q" — sq. Thus 

2 

g(s-)-oo) ~ Ai(s) - exp --s 3/2 (D.25) 

o 

It can be proved (lij that in the opposite limit, s — > — oo, the asymptotic form of the solution of the Panlevee equation 
(|D.24j) (which has the right tail Airy function limit, Eq. (|D.25jl ) is 



- x i ~ \/-i S (D.26) 



The Tracy- Widom distribution function Prwifi) is defined by the function ^(s) as follows. The function ^(s) 
gives the probability that a random quantity t described by a probability distribution functions Ptw (f ) bas the value 
less than a given parameter s: 



F 2 (s) 



dtP TW (t) 



(D.27) 



Taking the derivative of this relation and substituting here the result, Eq. (|D.4p , we find 



Ptw(s) = exp 



dt(t- s)q 2 {t) 



dtq 2 {t) 



(D.28) 



where the function q(s) is the solution of the differential equation (ID.24|) 
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Substituting the above two asymptotics into Eq. (|D.28p . we can estimate the asymptotic behavior for the right and 
the left tails of the TW probability distribution function: 



Ptw(s — > +00) ~ exp 
Ptw{ s —00) ~ exp 



-3 S 



12 



3/2 



(D.29) 
(D.30) 
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